El viento derivado de satélites se codifica como tipo 240 - 260 dependiendo del instrumento. Sin embargo ni GFS ni el resto de los sistemas operativos usan las observaciones presentes en el prepbufr y en cambio leen el bufr específico gdas.tHHz.satwnd.tm00.bufr_D.
En la Tabla 18 se detalla el uso de las distintas observaciones según el instrumento.
Mejor y con más detalle Tabla b
read_obs.f90
read_satwnd.f90
For the satellite ID type:
For satellite subtype:
The quality mark:the values range from 0 to 15.
Notas:
Prueba de asimilación usando 3DVar para hacer un par de experimentos rápidos.
Fecha: 2018112100
exp[bottom_top %in% seq(1, 60, 10)] %>%
ggplot(aes(west_east, south_north)) +
geom_point(aes(color = u)) +
scale_color_divergent() +
facet_grid(exp~bottom_top)
diff[bottom_top %in% seq(20, 58, 2)] %>%
ggplot(aes(west_east, south_north)) +
geom_point(aes(color = u.diff)) +
scale_color_divergent() +
facet_wrap(~bottom_top)
diag <- fread("/home/paola.corrales/comGSIv3.7_EnKFv1.3/examples/test/gsi_satwnd/results_conv_ges.2018112100") %>%
.[, c("V2", "V4") := NULL]
colnames(diag) <- c("var", "stationID", "type", "dhr", "lat", "lon", "pressure", "usage.flag", "obs", "obs.guess", "obs2", "obs.guess2", "rerr")
diag <- diag[type %between% c(240, 260)] %>%
.[, sat_type := case_when(
type == 240 ~ "GOES SW winds",
type == 241 ~ "India",
type == 242 ~ "JMA Visible",
type == 243 ~ "EUMETSAT visible",
type == 244 ~ "AVHRR winds",
type == 245 ~ "GOES IR",
type == 246 ~ "GOES WV cloud top",
type == 247 ~ "GOES WV deep layer",
type == 248 ~ "GOES cloud top (sounder)",
type == 249 ~ "GOES deep layer (sounder)",
type == 250 ~ "JMA WV deep layer",
type == 251 ~ "GOES visible",
type == 252 ~ "JMA IR winds",
type == 253 ~ "EUMETSAT IR winds",
type == 254 ~ "EUMETSAT WV deep layer winds",
type == 257 ~ "MODIS IR",
type == 258 ~ "MODIS WV cloud top",
type == 259 ~ "MODIS WV deep layer winds",
type == 260 ~ "VIIR IR winds")] %>%
melt(measure.vars = c("obs", "obs2", "obs.guess", "obs.guess2")) %>%
.[, var := if_else(str_detect(variable, "2"), "v", "u")] %>%
.[, variable := str_remove(variable, "2")] %>%
pivot_wider(names_from = variable, values_from = value) %>%
setDT
diag[var == "u"] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(aes(color = sat_type)) +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_wrap(~usage.flag, labeller = labeller(usage.flag = c("-1" = "rejected",
"1" = "asimilated"))) +
labs(color = "Source",
x = "lon") +
theme_minimal()
Comparando la distribución espacial de las observaciones rechazadas y asimiladas, no se ve gran diferencia. Lo único más notable es que todas las observaciones de EUMESAT son rechazadas, no queda claro si esto se debe a una casualidad (justo en este tipo estas observaciones eran de baja calidad o estaban muy lejos del guess) o si es algo definido en el código de GSI.
Revisando el código (read_satwnd.f90) parece que los “WV deep layer” se usan para monitoreo. También dice: “reject data zenith angle >68.0 degree”, habría que ver el ángulo zenital de estas observaciones. A simple vista hay muchas observaciones con un ángulo zenital menor a 68 grados así que esta no sería la causa para recharzar todas las observaciones.
Es notorio que el error de las observaciones rechazadas en la mayoría de los casos es Inf.
N_obs <- diag[var == "u", .N, by = .(usage.flag, var)] %>%
.[, label := paste("N = ", N)] %>%
.[, ":="(lon = 302.0,
pressure = 920)]
diag[var == "u"] %>%
ggplot(aes(ConvertLongitude(lon), pressure)) +
geom_point(aes(color = sat_type)) +
scale_y_level(name = "pressure level", breaks = c(1000, 900, 700, 500, 400, 300, 200, 100)) +
scale_x_longitude(breaks = seq(-75, -55, 5)) +
geom_label(data = N_obs, aes(x = ConvertLongitude(lon), y = pressure, label = label)) +
facet_grid(var~usage.flag, labeller = labeller(usage.flag = c("-1" = "rejected",
"1" = "asimilated"))) +
labs(color = "Source",
x = "lon",
y = "pressure level") +
theme_minimal()
diag %>%
ggplot((aes(ConvertLongitude(lon), lat))) +
geom_point(aes(color = obs)) +
scale_color_divergent() +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_grid(var~usage.flag, labeller = labeller(usage.flag = c("-1" = "rejected",
"1" = "asimilated"))) +
theme_minimal()
diag[var == "u"] %>%
ggplot(aes(ConvertLongitude(lon), pressure)) +
geom_point(aes(color = obs.guess)) +
scale_color_divergent() +
scale_y_level(name = "pressure level", breaks = c(1000, 900, 700, 500, 400, 300, 200, 100)) +
scale_x_longitude(breaks = seq(-75, -55, 5)) +
geom_label(data = N_obs, aes(x = ConvertLongitude(lon), y = pressure, label = label)) +
facet_grid(var~usage.flag, labeller = labeller(usage.flag = c("-1" = "rejected",
"1" = "asimilated"))) +
labs(color = "Source",
x = "lon",
y = "pressure level") +
theme_minimal()
10 miembros y 24 ciclos de asimilación arrancando a las 18 UTC del 20181120
files <- Sys.glob("analisis/diagfiles/E1/asim_conv*")
obs <- purrr::map(files, function(f) {
diag <- fread(f) %>%
.[, exp := basename(dirname(f))] %>%
.[, date := ymd_hms(stringr::str_extract(f, "\\d{14}"))] %>%
.[]
}) %>%
rbindlist() %>%
.[, c("V2", "V4") := NULL]
colnames(obs) <- c("var", "stationID", "type", "dhr", "lat", "lon", "pressure", "usage.flag", "flag.prep", "obs", "obs.guess", "obs2", "obs.guess2", "rerr", "exp", "date")
obs[, .N, by = .(exp, var, usage.flag)] %>%
dcast(exp + usage.flag ~ var, value.var = "N") %>%
knitr::kable()
| exp | usage.flag | ps | q | t | uv |
|---|---|---|---|---|---|
| E1 | -1 | 849979 | 4378 | 205425 | 1078653 |
| E1 | 1 | 171184 | 240339 | 939862 | 728592 |
obs[, .N, by = .(var, date, exp, usage.flag)] %>%
ggplot(aes(date, N)) +
geom_line(aes(color = var, linetype = factor(usage.flag))) +
scale_linetype_manual(values = c("-1" = 2, "1" = 1)) +
labs(title = "Cantidad de observaciones asimiladas y rechazadas según variable",
linetype = "Uso",
color = "experimento") +
theme_minimal()
Por defecto las observaciones tipo 251 = GOES Visible no se asimilan porque no tienen error definido. Algunas otras observaciones tienen defidino un error para algunos niveles y otros no. Esto permite controlar a que niveles se asimilacada observación, lo complicado sería redefinir estos errores para asimilar más o menos oservaciones sin conocer mejor el tipo de observación.
Respecto a la tabla convinfo, algunos tipos de observación tienen subtipos (algunos se asimilan, otros no) y no logro encontrar que significan canda uno. El gross check para estas observaciones es muy estricto con un umbral de entre 1.3 y 2.5 m/s (para comparar el umbral de las observaciones de estaciones de superficie es 6 m/s).
obs[var == "uv", .N, by = .(var, type, usage.flag)] %>%
dcast(usage.flag ~type, value.bar = "N") %>%
knitr::kable()
## Using 'N' as value column. Use 'value.var' to override
| usage.flag | 220 | 221 | 230 | 231 | 233 | 240 | 243 | 245 | 246 | 247 | 251 | 253 | 254 | 280 | 281 | 284 | 287 | 290 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| -1 | 4441 | 69 | NA | 97567 | 74 | 4334 | 2233 | 79544 | 71607 | 102012 | 203095 | 2863 | 10805 | NA | 17770 | 110 | 480740 | 1389 |
| 1 | 12664 | 371 | 22 | 16107 | 1873 | 737 | 242 | 42567 | 25215 | 83712 | 3188 | 206 | 327 | 77 | 223515 | NA | 298731 | 19038 |
obs <- obs[type %between% c(240, 260)] %>%
.[, sat_type := case_when(
type == 240 ~ "GOES SW winds",
type == 241 ~ "India",
type == 242 ~ "JMA Visible",
type == 243 ~ "EUMETSAT visible",
type == 244 ~ "AVHRR winds",
type == 245 ~ "GOES IR",
type == 246 ~ "GOES WV cloud top",
type == 247 ~ "GOES WV deep layer",
type == 248 ~ "GOES cloud top (sounder)",
type == 249 ~ "GOES deep layer (sounder)",
type == 250 ~ "JMA WV deep layer",
type == 251 ~ "GOES visible",
type == 252 ~ "JMA IR winds",
type == 253 ~ "EUMETSAT IR winds",
type == 254 ~ "EUMETSAT WV deep layer winds",
type == 257 ~ "MODIS IR",
type == 258 ~ "MODIS WV cloud top",
type == 259 ~ "MODIS WV deep layer winds",
type == 260 ~ "VIIR IR winds")]
obs[type %between% c(240, 260), .N, by = .(var,type, sat_type, date, exp, usage.flag)] %>%
ggplot(aes(date, N)) +
geom_line(aes(color = factor(type))) +
geom_point(aes(color = factor(type))) +
facet_wrap(usage.flag ~ exp) +
labs(title = "Observaciones de viento",
linetype = "Uso",
color = "type") +
theme_minimal()
Intentando idenficar la causa de rechazo de las observaciones veo que:
Infobs[, inf := !is.finite(rerr)] %>%
.[, .N, by = .(type, usage.flag, inf)] %>%
dcast(usage.flag + inf ~ type) %>%
knitr::kable()
## Using 'N' as value column. Use 'value.var' to override
| usage.flag | inf | 240 | 243 | 245 | 246 | 247 | 251 | 253 | 254 |
|---|---|---|---|---|---|---|---|---|---|
| -1 | FALSE | 1408 | 2098 | 57566 | 63620 | 66468 | 62648 | 339 | 2538 |
| -1 | TRUE | 2926 | 135 | 21978 | 7987 | 35544 | 140447 | 2524 | 8267 |
| 1 | FALSE | 737 | 242 | 42567 | 25215 | 83712 | 3188 | 206 | 327 |
source("help_functions.R")
path <- "analisis/diagfiles/E1/asim_conv_2018112*"
perfiles <- read_diag(path, variable = c("uv")) %>%
.[, c("error", "nivel.error") := input_obs_error(var, type, pressure)]
RCRV <- get_RCRV(perfiles[type %between% c(240, 260)], tipo = c("perfil")) %>%
.[, sat_type := case_when(
type == 240 ~ "GOES SW winds",
type == 241 ~ "India",
type == 242 ~ "JMA Visible",
type == 243 ~ "EUMETSAT visible",
type == 244 ~ "AVHRR winds",
type == 245 ~ "GOES IR",
type == 246 ~ "GOES WV cloud top",
type == 247 ~ "GOES WV deep layer",
type == 248 ~ "GOES cloud top (sounder)",
type == 249 ~ "GOES deep layer (sounder)",
type == 250 ~ "JMA WV deep layer",
type == 251 ~ "GOES visible",
type == 252 ~ "JMA IR winds",
type == 253 ~ "EUMETSAT IR winds",
type == 254 ~ "EUMETSAT WV deep layer winds",
type == 257 ~ "MODIS IR",
type == 258 ~ "MODIS WV cloud top",
type == 259 ~ "MODIS WV deep layer winds",
type == 260 ~ "VIIR IR winds")]
RCRV %>%
melt(measure.vars = c("mean.y", "sd.y")) %>%
ggplot(aes(nivel.error, value)) +
geom_hline(yintercept = c(0, 1), color = "darkgrey") +
geom_line(aes(color = factor(sat_type), linetype = variable)) +
scale_color_brewer(palette = "Set1") +
scale_x_level() +
coord_flip() +
facet_wrap(~ var, scales = "free_x") +
labs(title = "Reduced Centered Random Variable",
subtitle = "Satellite wind",
linetype = "",
color = "Obs type") +
theme_minimal()